В этом домашнем задании мы будем работать с датасетом Metabolic Syndrome, содержащим клинические и биохимические показатели пациентов.
Метаболический синдром - это комплекс взаимосвязанных нарушений обмена веществ, включающий ожирение, инсулинорезистентность, дислипидемию и артериальную гипертензию. Этот синдром значительно повышает риск развития сердечно-сосудистых заболеваний и диабета 2 типа.
Датасет включает следующие показатели:
Цель работы - применить методы многомерного анализа (PCA) и кластеризации для выявления структуры данных и классификации пациентов по наличию метаболического синдрома.
Зарегистрируйтесь на платформе Kaggle и скачайте датасет: Metabolic Syndrome Dataset
file <- 'data/Metabolic Syndrome.csv' # Обозначьте актуальный путь к файлу
# Загрузка данных
datama <- fread(file) %>%
select(-seqn, -Sex, -Marital, -Race, -Albuminuria, -Income)
# Удаление пропущенных значений
datama <- datama %>% na.omit()
# Сохранение целевой переменной
target <- datama$MetabolicSyndrome %>% as.factor() %>% recode('0' = 'Нет синдрома', '1' = 'Есть синдром')
datama$MetabolicSyndrome <- NULL
# Проверка размерности
cat("Размер датасета:", nrow(datama), "наблюдений,", ncol(datama), "переменных\n\n")## Размер датасета: 2311 наблюдений, 8 переменных
## Распределение классов:
## target
## Нет синдрома Есть синдром
## 1508 803
1 - наличие метаболитного синдрома 0 - отсутсвие
Для комплексного анализа взаимосвязей между всеми переменными
используем функцию ggpairs(), которая строит матрицу
диаграмм рассеяния, гистограмм и коэффициентов корреляции.
ggpairs(datama,
progress = F,
columnLabels = labels_name,
title = 'Матрица взаимосвязи переменных в исследовании метаболитного синдрома',
lower = list(continuous = wrap("points", colour = 'cyan3', alpha =0.3 ) )
)+
theme(
plot.title = element_text(size = 18),
axis.text.x = element_text(angle = 45)
)Подавляющее число переменных имеют статистически значимую, но слабую взаимосвязь (меньше 0,5). Статистически значимая (на уровне р<0,001) взаимосвязь обнаружена между ИМТ и окружностью талии (коэф. корр. Пирсона = 0,9). Статистически значимой взаимосвзи не обнаружено между ИМТ и уровнем альбумина в моче.
Вопрос: Какие переменные имеют правостороннюю асимметрию (скошены вправо, имеют длинный хвост вправо)?
Ответ: Правостороннюю ассиметрию имеют следующие пременные: уровень альбумина в моче, уровень глюкозы в крови, уровень триглециридов в крови, ИМТ, уровень холестерина.
Для переменных с сильной правосторонней асимметрией применяем логарифмическое преобразование, которое помогает: - Нормализовать распределение - Уменьшить влияние выбросов - Стабилизировать дисперсию
Применяем log1p() (логарифм с добавлением 1), чтобы
избежать проблем с нулевыми значениями:
Для корректной работы методов многомерного анализа необходимо стандартизировать данные (масштабирование к нулевому среднему и единичной дисперсии).
Вопрос: Какие переменные имеют аутлайеры (выбросы > 3SD или < -3SD)? Отразите диапазон (range) каждой переменной в шкалированном датасете в виде таблицы/матрицы.
dataoutliers <- data_scale %>%
as.data.frame() %>%
select(where(function(x) any(x>3 | x < -3)))
dataoutliers%>%
colnames()## [1] "WaistCirc" "BMI" "UrAlbCr" "UricAcid"
## [5] "BloodGlucose" "HDL" "Triglycerides"
dataoutliers %>%
pivot_longer(everything(), names_to = 'Variable', values_to = 'range') %>%
group_by(Variable) %>%
summarise(
Min = min(range) %>% round(2),
Max = max(range)%>% round(2)
) ## # A tibble: 7 × 3
## Variable Min Max
## <chr> <dbl> <dbl>
## 1 BMI -2.33 6.09
## 2 BloodGlucose -2.03 8.08
## 3 HDL -2.62 6.44
## 4 Triglycerides -1.07 15.1
## 5 UrAlbCr -1.4 5.69
## 6 UricAcid -2.57 4.07
## 7 WaistCirc -2.59 4.78
Ответ: Переменные с аутлайерами и их диапазоны значений: окружность талии (-2.59;4.78), ИМТ (-2.33;6.09), уровень альбумина в моче (-0,17;19,49), уровень мочевой кислоты в крови (-2.57;4.07), уровень глюкозы в крови (-2.03; 8.08), уровень холестерина (-2.62; 6.44), уровень триглециридов в крови (-1.07; 15.12).
Для уменьшения влияния экстремальных значений примените замену значений >3SD на 3 и <-3SD на -3:
ggpairs(data_scale_out,
progress = F,
columnLabels = labels_name,
title = 'Матрица взаимосвязи страндартизированных переменных в исследовании метаболитного синдрома',
lower = list(continuous = wrap("points", colour = 'cyan3', alpha =0.3 ) )) +
theme(
plot.title = element_text(size = 18)
)Метод главных компонент (PCA) - это линейное преобразование данных в новую систему координат, где: - Первая главная компонента объясняет максимальную долю вариации - Каждая последующая компонента объясняет максимум оставшейся вариации - Все компоненты ортогональны (некоррелированы) друг другу
## Standard deviations (1, .., p=8):
## [1] 1.5778616 1.1061183 0.9609338 0.8762480 0.7914170 0.6334983 0.6126126
## [8] 0.2625776
##
## Rotation (n x k) = (8 x 8):
## PC1 PC2 PC3 PC4 PC5
## Age 0.1989308 0.73901791 -0.01226988 0.09434503 0.55008085
## WaistCirc 0.5663908 -0.05530050 -0.38554347 -0.02372416 -0.01638634
## BMI 0.5092555 -0.14906467 -0.48402487 -0.06152298 -0.12179182
## UrAlbCr 0.1310311 0.49171471 0.20945961 -0.25023628 -0.76990209
## UricAcid 0.3525492 -0.01642545 0.36052040 0.82667532 -0.16508109
## BloodGlucose 0.2270847 0.21277497 0.14166833 -0.32603524 0.10311120
## HDL -0.3542312 0.37066253 -0.50154476 0.30535616 -0.16759452
## Triglycerides 0.2485549 -0.06210521 0.41407761 -0.20297648 0.15358535
## PC6 PC7 PC8
## Age -0.27246695 0.12797730 0.109645677
## WaistCirc -0.01544916 0.07526600 -0.721640170
## BMI 0.02140645 0.05481605 0.679762573
## UrAlbCr -0.17005463 0.11338931 -0.005569183
## UricAcid 0.09443564 -0.15488245 0.043746545
## BloodGlucose 0.49693855 -0.72050863 0.008828523
## HDL 0.58544914 0.14030725 -0.042606006
## Triglycerides 0.54546270 0.63189201 0.036029782
Биплот одновременно показывает: - Наблюдения (пациенты) - точки, окрашенные по диагнозу - Переменные - стрелки, указывающие направление максимального изменения
Постройте биплот используя fviz_pca_biplot(). Окрасьте точки по переменной target (диагноз)
pca_result_renam <- pca_result$var$coord
diagcol <- c('Нет синдрома' = 'mediumpurple',
'Есть синдром' = 'lightgoldenrod')
fviz_pca_biplot(pca_result,
geom = 'point',
habillage = 'none',
col.ind = target,
col.var = 'grey10',
repel = TRUE,
ggtheme = thememine
)+
scale_color_manual(
values = diagcol
)+
scale_x_continuous(breaks = seq(-4,5,1))+
scale_y_continuous(breaks = seq(-2,3,1))+
labs(
title = 'Биплот главных компонент',
x = 'Главная компонента 1',
y = 'Главная компонента 2',
colour = 'Наличие метаболитного синдрома',
shape = 'Наличие метаболитного синдрома'
)Вопрос 3.1: Какие две переменные вносят наибольший вклад в PC2?
Ответ: Наибольший вклад в PC2 вносят переменные Возраста и уровень альбумина, т.к. стрелки, обозначающие эти переменные имеют наибольшую длину проекции на вертикальную ось. Обе переменные положительно коррелируют с PC2.
Вопрос 3.2: Какие три переменные вносят наибольший вклад в PC1?
Ответ: наибольший вклад в PC1 вносят переменные ИМТ, измерение окружности тела, и почти на равном уровне вносят вклад переменные уровня мочевой кислоты в крови (положительная корреляция) и уровень холестерина (отрицательная корреляция). Если интерпертировать также матрицу нагрузок видно, что влияние уровня холестерина выше влияние уровня мочевой кислоты.
Вопрос 3.3: Опишите, как в пространстве главных компонент разделяются наблюдения на основе диагноза. Какие 5 переменных вносят наибольший вклад в разделение пациентов с метаболическим синдромом и без него?
Ответ: Переменные хорошо группируются по диагнозу вдоль горизонтальной оси, в обалсти положительных значений - пациенты с метаболическим синдромом, в области отрицательных - без. Наибольший вклад в PC1 и соответственно в разделение пациентов с метаболическим синдромом вносят следующие переменные: ИМТ, измерение окружности тела, уровень холестерина (отрицательная корреляция), уровень мочевой кислоты и уровень триглециридов в крови.
Вопрос 4.1: Сколько главных компонент необходимо для объяснения 90% вариации исходного датасета?
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.5779 1.1061 0.9609 0.8762 0.79142 0.63350 0.61261
## Proportion of Variance 0.3621 0.1779 0.1343 0.1117 0.09109 0.05836 0.05458
## Cumulative Proportion 0.3621 0.5400 0.6743 0.7859 0.87703 0.93539 0.98997
## PC8
## Standard deviation 0.26258
## Proportion of Variance 0.01003
## Cumulative Proportion 1.00000
Ответ: Для объяснения не менее 90% вариаций исходного датасета необходиом 6 копомнент, которые объяснят 93,5% вариаций.
Вопрос 4.2: Какой процент вариации объясняют первые 3 главные компоненты?
Ответ: Первые три главные компоненты объяснят 67,4% вариаций.
Метод локтя помогает определить, после какой компоненты добавление новых компонент дает незначительный прирост объясненной вариации.
# Ваш код здесь
fviz_eig(pca_result,
addlabels = TRUE,
ncp = 8,
main = "Scree plot: доля объясненной дисперсии",
ggtheme = thememine)+
labs(
x = 'Главные компоненты',
y = 'Доля объясненной дисперсии'
)+
scale_y_continuous(breaks = seq(0,40,5))+
theme(
plot.title = element_text(size = 18),
axis.text = element_text(size = 12),
axis.title = element_text(size = 15)
)Вопрос: Используя метод локтя, определите оптимальное количество главных компонент. Какой процент вариации они объясняют?
Ответ: Оптимальное по методу локты количество главные компонент - 3, которые объясняют 64,7% вариации.
Метод Кайзера (Kaiser criterion) основан на собственных значениях (eigenvalues).
Kaiservalues = pca_result$sdev ^ 2
names(Kaiservalues) = paste0("PC", 1:length(Kaiservalues))
Kaiservalues %>% round(2)## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 2.49 1.22 0.92 0.77 0.63 0.40 0.38 0.07
Вопрос: Сколько главных компонент метод Кайзера рекомендует использовать для дальнейшего анализа?
Ответ: По методу Кайзера можно использовать 2 главные компоненты для дальнейшего анализа.
Практические рекомендации:
Для нашего анализа: Мы будем исследовать варианты с разным количеством компонент и выберем оптимальное на основе качества кластеризации.
Для понимания биологического смысла главных компонент необходимо проанализировать, какие исходные переменные вносят наибольший вклад в каждую компоненту.
Вычислите корреляции между исходными переменными и главными компонентами. Создайте тепловую карту используя pheatmap().
сorrlvar <- cor(data_scale_out, pca_result$x)
pheatmap(сorrlvar,
cluster_rows = TRUE,
cluster_cols = FALSE,
display_numbers = TRUE,
number_format = "%.2f",
main = "Корреляция переменных с главными компонентами"
)Вопрос: Какие две переменные наиболее сильно коррелируют с каждой из первых трех главных компонент (PC1, PC2, PC3)?
Ответ:
PC1:
- Переменная 1: [окружность тела] (корреляция: [0,90])
- Переменная 2: [ИМТ] (корреляция: [0,83])
PC2:
- Переменная 1: [возраст] (корреляция: [0,82])
- Переменная 2: [уровень альбумина в моче] (корреляция: [0,60])
PC3:
- Переменная 1: [уровень холестерина и уровень триглицеридов в крови] (корреляция: [-0,50 и 0,50 соответсвенно])
- Переменная 2: [ИМТ] (корреляция: [-0,48])
Примените иерархическую кластеризацию, постройте матрицу сопряженности и оцените совпадение предсказанных классов и диагноза с помощью Хи-квадрата.
####################### Параметры
numPCs <- 6
distMethod = 'euclidean' # #'manhattan'
hclustMethod = 'ward.D2' # 'average' # # #//'complete'
#######################
PCs = pca_result$x[,1:numPCs]
hclust_pc6 <-
PCs %>%
dist(method = distMethod) %>%
hclust(method = hclustMethod)
hclusters_pc6_2 <- cutree(hclust_pc6, 2)
# Посчитаем матрицу сопряженности
contingency_table <- table(Кластер = hclusters_pc6_2, Диагноз = target)
contingency_table## Диагноз
## Кластер Нет синдрома Есть синдром
## 1 861 757
## 2 647 46
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: contingency_table
## X-squared = 343.16, df = 1, p-value < 2.2e-16
Для построения кластеризации использовали 6 главные компонент, объясняющих 93,5% вариаций, евклидово расстояние как мера расчета расстояний, т.к. наши данных находятся в одной размерности, и ward.D как метод иерархической кластеризации, инимизирующий внутрикластерную дисперсию. при котором значения X-squared наивысшие среди трех остальных.
Согласно полученному критерию Хи-квадрат для проверки независимости, разделение на кластеры и поставленный диагноз имеют статистически значимую взаимосвязь по уровню значимости 0.05 (p < 2.2e-16, X-квадрат = 469.93).
При этом, кластер 1 на 87,5% состоит из пациентов, не имеющих метаболического синдрома, а кластер 2 на 70,0% из пациентов, имеющих метаболический синдром.
Протестируйте различные комбинации и выберите ту которая даёт лучшую кластеризацию на основне Хи-квадрата: - Количество главных компонент: 1-6 (numPCs) - Метрики расстояния: euclidean, manhattan (distMethod) - Методы связывания: ward.D2, complete, average (hclustMethod)
Статистика теста Хи-квадрата наибольшая при использовании одной компоненты. Остальные компоненты возможно вносят излишную информацию, не так связаную с постановлением диагноза. Изменение меры расчета расстояния не влияет на получаемую статистику. Оба метода равноэффективны, т.к. наши данные находся в одной размерности и не имеют аутлайнеров. ward.D как метод иерархической кластеризации дает наибольшие значения X-squared среди трех остальных методов.
Создайте тепловую карту с иерархической кластеризацией. Используйте оптимальные параметры. Добавьте аннотацию с диагнозами.
Тепловая карта с использованием всех главных компонент. Кластер с большинством пациентов с синдромом также имеет много наблюдений пациентов без синдрома.
# Создание фрейма данных с аннотациями
annotation_row <-
data.frame(
'Диагноз' = target)
row.names(data_scale_out) <- paste0('obs', 1:nrow(data_scale_out))
row.names(annotation_row) <- row.names(data_scale_out)
# Тепловая карта с аннотациями кластеров
pheatmap(data_scale_out,
# Кластеризация
cluster_cols = F,
cluster_rows = T,
clustering_method = 'ward.D2',
clustering_distance_rows = 'euclidean',
cutree_rows = 2,
# Внешний вид
show_rownames = F, # Убираем подписывание колонок (индексы строк/наблюдений),
angle_col = 45, # Подписи под углом
# Аннотации
annotation_row = annotation_row,
annotation_colors = list('Диагноз' = c('Нет синдрома' = "mediumpurple",
'Есть синдром' = "lightgoldenrod"))
)При построении тепловой карты с использованием одной компоненты, как оптимальным количеством (по результатам Хи-квадрат), визуально лучше разделяет пациентов. В кластере с преимущественным расположением пациентов без синдрома располагается больше наблюдений пациентов с синдромом, но второй кластер по сравнению с первым гарфиком преимущественно состоит из пациентов с синдромом.
# Создание фрейма данных с аннотациями
annotation_row <-
data.frame(
'Диагноз' = target)
row.names(data_scale_out) <- paste0('obs', 1:nrow(data_scale_out))
row.names(annotation_row) <- row.names(data_scale_out)
pca1mart <- data.frame(PC1 = pca_result$x[,1])
rownames(pca1mart) <- paste0('obs', 1:nrow(data_scale_out))
pca1mart <- as.matrix(pca1mart)
# Тепловая карта с аннотациями кластеров
pheatmap(pca1mart,
# Кластеризация
cluster_cols = F,
cluster_rows = T,
clustering_method = 'ward.D2',
clustering_distance_rows = 'euclidean',
cutree_rows = 2,
# Внешний вид
show_rownames = F, # Убираем подписывание колонок (индексы строк/наблюдений),
angle_col = 45, # Подписи под углом
# Аннотации
annotation_row = annotation_row,
annotation_colors = list('Диагноз' = c('Нет синдрома' = "mediumpurple",
'Есть синдром' = "lightgoldenrod"))
)Теперь сравним результаты с методом K-means, который использует другой подход к кластеризации.
Протестируйте K-means с разным количеством главных компонент (1-6).
# Ваш код здесь
PCs = pca_result$x[,1:2]
kmeans_pc6_2 = kmeans(PCs, centers = 2, nstart = 25)
contingency_table <- table(Кластер = kmeans_pc6_2$cluster, Диагноз = target)
contingency_table## Диагноз
## Кластер Нет синдрома Есть синдром
## 1 364 662
## 2 1144 141
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: contingency_table
## X-squared = 719.16, df = 1, p-value < 2.2e-16
Вопрос: Какое минимальное количество PCs даёт наилучшую кластеризацию основываясь на Хи-квадрат?
Ответ: 2 главные компоненты дают наибольшую статистику Хи-квадрат.
Вопрос: На основе статистики χ² (Chi-squared), какой метод и с каким количеством главных компонент дает наилучший результат классификации?
Ответ: Метод методом K-means с двумя главными компонетами дает значения Хи-квадрат несколько выше значений, полученных с помощью иерархической кластеризацией с использованием 1 компоненты.
Постройте два биплота fviz_pca_biplot() рядом:
1) С окраской по реальному диагнозу
2) С окраской по кластерам лучшего метода
Используйте gridExtra::grid.arrange(…, ncol = 2) для размещения
графиков
# Ваш код здесь
plottarg <- fviz_pca_biplot(pca_result,
geom = 'point',
habillage = 'none',
col.ind = target,
col.var = 'grey10',
repel = TRUE,
ggtheme = thememine
)+
scale_color_manual(
values = diagcol
)+
scale_x_continuous(breaks = seq(-4,5,1))+
scale_y_continuous(breaks = seq(-2,3,1))+
labs(
title = 'Биплот главных компонент',
x = 'Главная компонента 1',
y = 'Главная компонента 2',
colour = 'Наличие метаболитного синдрома',
shape = 'Наличие метаболитного синдрома'
)+
theme(
legend.position = 'bottom',
legend.justification = 'centre'
)
plotclus <- fviz_pca_biplot(pca_result,
geom = 'point',
habillage = 'none',
col.ind = factor(kmeans_pc6_2$cluster),
col.var = 'grey10',
repel = TRUE,
ggtheme = thememine
)+
scale_x_continuous(breaks = seq(-4,5,1))+
scale_y_continuous(breaks = seq(-2,3,1))+
scale_colour_manual(
values = c('1' = 'plum1',
'2' = 'bisque3')
)+
labs(
title = 'Биплот главных компонент',
x = 'Главная компонента 1',
y = 'Главная компонента 2',
colour = 'Кластер',
shape = 'Кластер'
)+
theme(
legend.position = 'bottom',
legend.justification = 'centre'
)+
theme(
legend.position = 'bottom',
legend.justification = 'centre'
)
gridExtra::grid.arrange(plottarg, plotclus, ncol = 2)На этом заканчивается основная часть задания! Выполнив её, вы сможете получить 100 баллов!
# Настройки UMAP
umap_config <- umap.defaults
umap_config$n_neighbors <- 50
umap_config$min_dist <- 0.01
umap_config$metric <- 'euclidean'
umap_config$random_state <- 42
# Применяем UMAP
umap_result <- umap(data_scale_out, config = umap_config)
# Визуализация
umap_data <- data.frame(
UMAP1 = umap_result$layout[, 1],
UMAP2 = umap_result$layout[, 2],
Diagnosis = target
)
ggplot(umap_data, aes(x = UMAP1, y = UMAP2, color = Diagnosis)) +
geom_point(alpha = 0.6, size = 2) +
scale_color_manual(
name = 'Наличие метаболитного синдрома',
values = diagcol) +
scale_x_continuous(breaks = seq(-4,4,1))+
scale_y_continuous(breaks = seq(-3,3,1))+
labs(title = "UMAP визуализация",
x = "UMAP 1",
y = "UMAP 2") +
theme() +
theme(legend.position = "right")Из графика видно, что группы частично разделяются, пациенты с синдромом сконцентрированы в верхней части графика (UMAP2>1), без синдрома (UMAP2<-1) - в нижней. Однко из-за схожести групп кластеризация не идеальная и группы перемешиваются между собой преимущесвенно в центре графика и частично у полюсов второго кластера.